dosresmeta R packageEvidence Synthesis and Meta-Analysis in R, March 27-31, 2023
A vignette of the dosresmeta package
Derive a dose-response curve from published results of multiple studies on the relation between a quantitative exposure (e.g. number of cigarettes, diet) and the occurrence of a health outcome (e.g. mortality risk, incidence of a disease).
Research questions:
Is there any association? What is the shape of the association?
What are the exposure values (doses) associated with the best or worst outcome?
What are the factors that can influence the dose–response shape?
Several fields of application.
Many leading medical and epidemiological journals.
International health organizations and academic institutions.
Measures of public health impact.
library(dosresmeta)
library(tidyverse)
data(coffee_mort)
gt(filter(coffee_mort, id < 3), groupname_col = c("id", "author"))| year | type | dose | cases | n | logrr | se | gender | area |
|---|---|---|---|---|---|---|---|---|
| 1 - LeGrady et al. | ||||||||
| 1987 | ci | 0.5 | 57 | 249 | 0.0000000 | 0.0000000 | M | USA |
| 1987 | ci | 2.5 | 136 | 655 | -0.2876821 | 0.1391187 | M | USA |
| 1987 | ci | 4.5 | 144 | 619 | -0.1743534 | 0.1373198 | M | USA |
| 1987 | ci | 6.5 | 115 | 387 | 0.0861777 | 0.1401409 | M | USA |
| 2 - Rosengren et al. | ||||||||
| 1991 | ci | 0.0 | 17 | 192 | 0.0000000 | 0.0000000 | M | Europe |
| 1991 | ci | 1.5 | 88 | 1121 | -0.1203576 | 0.2531537 | M | Europe |
| 1991 | ci | 3.5 | 155 | 2464 | -0.3418342 | 0.2442559 | M | Europe |
| 1991 | ci | 5.5 | 155 | 1986 | -0.1261707 | 0.2440559 | M | Europe |
| 1991 | ci | 7.5 | 39 | 575 | -0.2665264 | 0.2784189 | M | Europe |
| 1991 | ci | 9.5 | 24 | 427 | -0.5108256 | 0.2802582 | M | Europe |
p_1s <- subset(coffee_mort, id == 1) %>%
mutate(
low_rr = exp(logrr - qnorm(.975)*se),
upp_rr = exp(logrr + qnorm(.975)*se)) %>%
ggplot(aes(dose, exp(logrr))) +
geom_point() +
geom_errorbar(aes(ymin = low_rr, ymax = upp_rr),
col = "grey", width = .1) +
scale_y_continuous(trans = "log") +
labs(x = "Coffee consumption, cups/day",
y = "Relative risk",
title = "LeGrady et al., 1987")Observations are not independent (same reference category).
The RR in the referent is 1 by definition (log RR = 0).
Log-linear model: \(y = X \beta + \epsilon\)
\(y\) : vector of non-referent log \(RR\).
\(X\) : design matrix contains the assigned doses.
Additional info: number of cases, and number of subjects or person-time.
Assuming a linear trend: \(y = \beta_1 x + \epsilon\)
mods_lin <- dosresmeta(formula = logrr ~ dose, se = se, type = "ci",
cases = cases, n = n, data = subset(coffee_mort, id == 1))
Epi::ci.exp(mods_lin) exp(Est.) 2.5% 97.5%
dose 1.033376 0.9914259 1.077101
Every 1 cup/day of coffee increase is associated with a 3% (95% CI: 0.99, 1.07) higher mortality risk.
library(rms)
k <- quantile(coffee_mort$dose, c(.1, .5, .9))
mods_spl <- dosresmeta(formula = logrr ~ rcs(dose, k), se = se, type = "ci",
cases = cases, n = n, data = subset(coffee_mort, id == 1))
summary(mods_spl)Call: dosresmeta(formula = logrr ~ rcs(dose, k), type = "ci", cases = cases,
n = n, data = subset(coffee_mort, id == 1), se = se)
One-stage fixed-effects meta-analysis
Covariance approximation: Greenland & Longnecker
Chi2 model: X2 = 11.5464 (df = 2), p-value = 0.0031
Fixed-effects coefficients
Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
rcs(dose, k)dose -0.2048 0.0814 -2.5156 0.0119 -0.3644 -0.0452 *
rcs(dose, k)dose' 0.3930 0.1300 3.0225 0.0025 0.1382 0.6479 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
1 study 3 values, 2 fixed and 0 random-effects parameters
logLik AIC BIC
3.7764 -3.5527 -5.3555
xref <- .5
preds_spl <- data.frame(dose = c(xref, seq(0, 7, .1))) %>%
cbind(spl = predict(mods_spl, newdata = ., expo = T), lin = predict(mods_lin, newdata = ., expo = T))
p_1s + geom_line(data = preds_spl, aes(y = spl.pred, linetype = "Curve Spline")) +
geom_line(data = preds_spl, aes(y = lin.pred, linetype = "Linear")) +
labs(linetype = "Model")p_all <- coffee_mort %>% mutate(coffee_mort, low_rr = exp(logrr - qnorm(.975)*se), upp_rr = exp(logrr + qnorm(.975)*se)) %>%
ggplot(aes(dose, exp(logrr))) + geom_point() +
geom_errorbar(aes(ymin = low_rr, ymax = upp_rr), col = "grey", width = .1) +
facet_wrap(~ id, scales = "free") + scale_y_continuous(trans = "log", breaks = pretty_breaks()) +
labs(x = "Coffee consumption, cups/day", y = "Relative risk")
p_allApply the same dose-response analyses across studies.
And get model-based predictions.
predi <- lapply(unique(coffee_mort$id), function(i) {
xref <- coffee_mort$dose[coffee_mort$id == i & coffee_mort$se == 0]
data.frame(dose = c(xref, seq(0, max(coffee_mort$dose[coffee_mort$id == i]), .1))) %>%
cbind(spl = predict(modi[[as.character(i)]]$spl, newdata = ., expo = T),
lin = predict(modi[[as.character(i)]]$lin, newdata = ., expo = T))
}) %>% set_names(nm = unique(coffee_mort$id)) %>%
bind_rows(.id = "id")
p_all <- p_all + geom_line(data = predi, aes(y = spl.pred, linetype = "Curve Spline")) +
geom_line(data = predi, aes(y = lin.pred, linetype = "Linear")) +
labs(linetype = "Model")
p_allContrast and combine the dose-response coefficients \(b_i\).
\[b_i \sim N \left( \bar b , v_{b_i} + \tau^2 \right)\]
\[\hat {\bar b} = \frac{\sum_{i = 1}^K b_i w_i}{\sum_{i = 1}^K w_i}\] \(w_i = (v_{b_i} + \tau^2)^{-1}\)
\(\tau^2\) is the between-studies heterogeneity.
Test for overall association: \(H_0: \bar b = 0\).
Test (Cochran \(Q\)) and quantify the impact (\(I^2\)) of statistical heterogeneity.
id yi vi
1.dose 1 0.03283124 4.470869e-04
2.dose 2 -0.02360445 3.561441e-04
3.dose 3 -0.01430817 5.833618e-05
4.dose 4 -0.04777017 6.194806e-04
5.dose 5 -0.04736154 1.219069e-03
6.dose 6 -0.02027627 1.127937e-04
Random-Effects Model (k = 22; tau^2 estimator: REML)
logLik deviance AIC BIC AICc
37.5674 -75.1348 -71.1348 -69.0457 -70.4681
tau^2 (estimated amount of total heterogeneity): 0.0003 (SE = 0.0002)
tau (square root of estimated tau^2 value): 0.0173
I^2 (total heterogeneity / total variability): 75.01%
H^2 (total variability / sampling variability): 4.00
Test for Heterogeneity:
Q(df = 21) = 77.0088, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
-0.0326 0.0051 -6.4424 <.0001 -0.0425 -0.0227 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
id yi.1 yi.2 vi.1 vi.2 vi.3 vi.4
1 1 -0.20483274 0.39304403 0.0066300437 -0.0102252526 -0.0102252526 0.016910322
2 2 -0.07391787 0.07231451 0.0071693533 -0.0097924956 -0.0097924956 0.014074567
3 3 -0.02578942 0.02177087 0.0005222641 -0.0008797053 -0.0008797053 0.001668107
4 4 -0.17681455 0.30128718 0.0053709253 -0.0110934653 -0.0110934653 0.025900538
5 5 -0.08529809 0.08507452 0.0107967876 -0.0214784900 -0.0214784900 0.048166538
6 6 -0.20065255 0.24137081 0.0034649927 -0.0044857506 -0.0044857506 0.006002614
Call: mvmeta(formula = cbind(yi.1, yi.2) ~ 1, S = bi_spl[, c(4, 5,
7)], data = bi_spl)
Multivariate random-effects meta-analysis
Dimension: 2
Estimation method: REML
Fixed-effects coefficients
Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
yi.1 -0.0811 0.0115 -7.0780 0.0000 -0.1035 -0.0586 ***
yi.2 0.1072 0.0217 4.9451 0.0000 0.0647 0.1497 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Between-study random-effects (co)variance components
Structure: General positive-definite
Std. Dev Corr
yi.1 0.0372 yi.1
yi.2 0.0707 -0.9601
Multivariate Cochran Q-test for heterogeneity:
Q = 101.3912 (df = 42), p-value = 0.0000
I-square statistic = 58.6%
22 studies, 44 observations, 2 fixed and 3 random-effects parameters
logLik AIC BIC
46.2962 -82.5924 -73.9040
dosresmeta packagedosresmeta?It facilitates both the study-specific dose-response and meta-analysis.
It gives a comprehensive framework for dose-response meta-analysis.
It provides functions for testing, predicting, and interpreting (presenting) results.
It implements extensions (one-stage procedure, continuous outcome).
mod_lin <- dosresmeta(formula = logrr ~ dose, id = id, se = se, type = type,
cases = cases, n = n, data = coffee_mort, method = "mm")
summary(mod_lin)Call: dosresmeta(formula = logrr ~ dose, id = id, type = type, cases = cases,
n = n, data = coffee_mort, se = se, method = "mm")
Two-stage random-effects meta-analysis
Estimation method:
Covariance approximation: Greenland & Longnecker
Chi2 model: X2 = 44.4226 (df = 1), p-value = 0.0000
Fixed-effects coefficients
Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
(Intercept) -0.0324 0.0049 -6.6650 0.0000 -0.0419 -0.0229 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Between-study random-effects (co)variance components
Std. Dev
0.0163
Univariate Cochran Q-test for residual heterogeneity:
Q = 77.0088 (df = 21), p-value = 0.0000
I-square statistic = 72.7%
22 studies, 22 values, 1 fixed and 1 random-effects parameters
mod_spl <- dosresmeta(formula = logrr ~ rcs(dose, k), id = id, se = se, type = type,
cases = cases, n = n, data = coffee_mort, method = "mm")
summary(mod_spl)Call: dosresmeta(formula = logrr ~ rcs(dose, k), id = id, type = type,
cases = cases, n = n, data = coffee_mort, se = se, method = "mm")
Two-stage random-effects meta-analysis
Estimation method:
Covariance approximation: Greenland & Longnecker
Chi2 model: X2 = 209.5069 (df = 2), p-value = 0.0000
Fixed-effects coefficients
Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
rcs(dose, k)dose -0.0764 0.0105 -7.3101 0.0000 -0.0969 -0.0559 ***
rcs(dose, k)dose' 0.1053 0.0229 4.5935 0.0000 0.0604 0.1502 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Between-study random-effects (co)variance components
Std. Dev Corr
rcs(dose, k)dose 0.0354 rcs(dose, k)dose
rcs(dose, k)dose' 0.0797 -1
Univariate Cochran Q-test for residual heterogeneity:
Q = 101.3912 (df = 42), p-value = 0.0000
I-square statistic = 58.6%
22 studies, 44 values, 2 fixed and 1 random-effects parameters
data.frame(dose = c(0, seq(0, 6, .1))) %>%
cbind(spl = predict(mod_spl, newdata = ., expo = T), lin = predict(mod_lin, newdata = ., expo = T)) %>%
ggplot(aes(dose, spl.pred)) +
geom_line(aes(linetype = "Curve Spline")) +
geom_ribbon(aes(ymin = spl.ci.lb, ymax = spl.ci.ub), alpha = .1) +
scale_y_continuous(trans = "log") +
labs(x = "Coffee consumption, cups/day", y = "Relative risk", linetype = "Model")The normality assumption on the random-effects can be exploited to better predict study-specific curves.
The empirical Bayes estimates’ or ‘Best Linear Unbiased Predictions’ can be computed as:
\[\hat b_i = \hat {\bar b} + \tau^2 v_{b_i}^{-1}(b_i - \hat {\bar b})\]
# study-specific dose-response coefficients (blup)
bi_blup <- t(apply(blup(mod_spl), 1, function(x) x + rbind(coef(mod_spl)))) %>%
cbind(id = as.double(levels(mod_spl$id)))
predi_blup <- lapply(unique(coffee_mort$id), function(i) {
xref <- coffee_mort$dose[coffee_mort$id == i & coffee_mort$se == 0]
newdata <- data.frame(dose = c(xref, seq(0, max(coffee_mort$dose[coffee_mort$id == i]), .1)))
ctr.mat <- rcspline.eval(newdata$dose, k, inclx = T)
ctr.mat <- scale(ctr.mat, center = ctr.mat[1, ], scale = FALSE) # important to center the variables at the ref
data.frame(dose = newdata$dose, blup = exp(ctr.mat %*% bi_blup[bi_blup[, "id"] == i, 1:2]))
}) %>% set_names(nm = unique(coffee_mort$id)) %>%
bind_rows(.id = "id")
p_all + geom_line(data = predi_blup, aes(y = blup, linetype = "BLUP"), col = "blue")Goodness-of-fit statistics:
Deviance test:
D = 141.332 (df = 77), p-value = 0.000
Coefficient of determination R-squared: 0.679
Adjusted R-squared: 0.671
To further explain heterogeneity and/or explore differences in dose-response curves.
Prespecified study-level characteristics.
Limited number of data points: interpret with caution.
Stratified (or subgroups) analysis: divide the observations in groups and evaluate differences.
Meta-regression: formalization of stratified analysis, + investigation of continuous variables (be aware of pitfall).
spl_regr <- dosresmeta(formula = logrr ~ rcs(dose, k), id = id, se = se, type = type,
cases = cases, n = n, data = coffee_mort, mod = ~ area)
spl_regrCall: dosresmeta(formula = logrr ~ rcs(dose, k), id = id, type = type,
cases = cases, n = n, data = coffee_mort, mod = ~area, se = se)
Fixed-effects coefficients:
rcs(dose, k)dose.(Intercept) rcs(dose, k)dose'.(Intercept) rcs(dose, k)dose.areaJapan rcs(dose, k)dose'.areaJapan rcs(dose, k)dose.areaUSA rcs(dose, k)dose'.areaUSA
-0.0898 0.1014 -0.0497 0.1573 0.0221 -0.0073
22 studies, 44 values, 6 fixed and 3 random-effects parameters
logLik AIC BIC
45.9634 -73.9267 -59.1884
Is there any evidence of differential association across geographical areas?
expand.grid(dose = seq(0, 7, .2), area = levels(coffee_mort$area)) %>%
cbind(predict(spl_regr, newdata = ., expo = T)) %>%
ggplot(aes(x = dose, y = pred, col = area)) +
geom_line() +
scale_y_continuous(trans = "log") +
labs(x = "Coffee consumption, cups/day", y = "Relative risk", linetype = "Model")Unified approach which combines dose-response and pooling estimation (linear mixed models for summarized data).
Advantage:
# reconstructing covariances
Slist <- split(coffee_mort_all, coffee_mort_all$id) %>% lapply(function(x) {
if (any(is.na(x$cases) | is.na(x$n))) {
diag(x$se[x$se != 0 & !is.na(x$se)]^2, nrow = sum(x$se != 0 & !is.na(x$se)))
} else {
covar.logrr(y = logrr, v = I(se^2), cases = cases, n = n,
type = type, data = x)
}
})
# Splines without exclusion
spl_noexcl <- dosresmeta(logrr ~ rcs(dose, k), id = id, type = type,
cases = cases, n = n, se = se, data = coffee_mort_all,
covariance = "user", Slist = Slist, proc = "1stage")# Graphical comparison
data.frame(dose = c(0, seq(0, 6, .1))) %>%
cbind(spl_2s = predict(mod_spl, newdata = ., expo = T),
spl_1s = predict(spl_noexcl, newdata = ., expo = T)) %>%
ggplot(aes(dose)) +
geom_line(aes(y = spl_2s.pred, linetype = "Two-stage")) +
geom_line(aes(y = spl_1s.pred, linetype = "One-stage")) +
scale_y_continuous(trans = "log") +
labs(x = "Coffee consumption, cups/day", y = "Relative risk",
linetype = "Approach")Extension of the methodology where the outcome is expressed as means/differences in means instead of (log) relative risks.
https://alecri.github.io/software/dosresmeta.html Reproducible code, Analysis examples, and Additional useful codes.
https://github.com/alecri with more reproducible codes.